{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "![Curso Schwarz-Sosa-Suriano](http://www.fi.uba.ar/sites/default/files/logo.png)\n", "# Integración Numérica - Tercera parte\n", "**Curso Schwarz - Sosa - Suriano**\n", "- Métodos Numéricos. *Curso 2*\n", "- Análisis Numérico I. *Curso 4*\n", "- Métodos Matemáticos y Numéricos. *Curso 6*" ] }, { "cell_type": "markdown", "metadata": { "hideCode": true, "slideshow": { "slide_type": "slide" } }, "source": [ "## Cuadratura de Gauss-Legendre\n", "\n", "Este método propone una perspectiva totalmente distinta a la de los métodos de Newton-Cotes. En lugar de simplificar la elección de los $x_i$ tomando puntos equidistantes para obtener un sistema de ecuaciones lineales y obtener fácilmente los los $c_i$, *mantiene la totalidad de las incógnitas como tales con la finalidad de obtener una mejora en el orden del error*." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "Para generar el Sistema de Ecuaciones Lineales de $2.n+2$ incógnitas y $2.n+2$ incógnitas, propone primero cambiar la variable $x \\in [a,b]$ por $t \\in [-1,1]$:\n", "\n", "$$ \\int_a^b f(x) = \\frac{b-a}{2} \\int_{-1}^1 f\\left(\\frac{b-a}{2}t \n", "+ \\frac{a+b}{2}\\right)\\,dt$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Así las cosas, nuestra fórmula de cuadratura queda como:\n", "\n", "$$ Q_n = \\frac{b-a}{2} \\sum_{i=0}^n c_i . f\\left(\\frac{b-a}{2}t_i + \\frac{a+b}{2}\\right)$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Esta vez, al buscar las $2.n+2$ ecuaciones se aseguran resultados exactos al integrar polinomios de hasta grado $2.n+1$ *-lo cual implica además una mejora en el término de error de truncamiento.*" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "De la resolución del Sistema de Ecuaciones No Lineales para este nuevo intervalo normalizado surge que los $t_i$ coinciden con las raíces de los polinomios de Legendre, y los valores de los coeficientes $c_i$ se encuentran también relacionados con los mismos." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Aunque es posible calcular ambos valores en forma analítica, para aplicar el método se suele recurrir a los valores tabulados precalculados. De hecho, los tenemos disponibles en Numpy:" ] }, { "cell_type": "code", "execution_count": 455, "metadata": { "scrolled": false, "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "tTabla, cTabla, pTabla, nTabla = [], [], [], []\n", "for i in range (1,8):\n", " x, c = np.polynomial.legendre.leggauss(i)\n", " ok = False\n", " for k in range (len(x)):\n", " tTabla.append(x[k])\n", " cTabla.append(c[k])\n", " if (not ok):\n", " nTabla.append(i-1)\n", " pTabla.append(str(i)+' Puntos')\n", " ok = True\n", " else:\n", " pTabla.append('')\n", " nTabla.append('')\n", "\n", "df = pd.DataFrame(data = {'$n$': nTabla, '$t_i$': tTabla, '$c_i$': cTabla}, index= pTabla) " ] }, { "cell_type": "code", "execution_count": 457, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
$n$$t_i$$c_i$
1 Puntos00.0000002.000000
2 Puntos1-0.5773501.000000
0.5773501.000000
3 Puntos2-0.7745970.555556
0.0000000.888889
0.7745970.555556
4 Puntos3-0.8611360.347855
-0.3399810.652145
0.3399810.652145
0.8611360.347855
5 Puntos4-0.9061800.236927
-0.5384690.478629
0.0000000.568889
0.5384690.478629
0.9061800.236927
6 Puntos5-0.9324700.171324
-0.6612090.360762
-0.2386190.467914
0.2386190.467914
0.6612090.360762
0.9324700.171324
7 Puntos6-0.9491080.129485
-0.7415310.279705
-0.4058450.381830
0.0000000.417959
0.4058450.381830
0.7415310.279705
0.9491080.129485
\n", "
" ], "text/plain": [ " $n$ $t_i$ $c_i$\n", "1 Puntos 0 0.000000 2.000000\n", "2 Puntos 1 -0.577350 1.000000\n", " 0.577350 1.000000\n", "3 Puntos 2 -0.774597 0.555556\n", " 0.000000 0.888889\n", " 0.774597 0.555556\n", "4 Puntos 3 -0.861136 0.347855\n", " -0.339981 0.652145\n", " 0.339981 0.652145\n", " 0.861136 0.347855\n", "5 Puntos 4 -0.906180 0.236927\n", " -0.538469 0.478629\n", " 0.000000 0.568889\n", " 0.538469 0.478629\n", " 0.906180 0.236927\n", "6 Puntos 5 -0.932470 0.171324\n", " -0.661209 0.360762\n", " -0.238619 0.467914\n", " 0.238619 0.467914\n", " 0.661209 0.360762\n", " 0.932470 0.171324\n", "7 Puntos 6 -0.949108 0.129485\n", " -0.741531 0.279705\n", " -0.405845 0.381830\n", " 0.000000 0.417959\n", " 0.405845 0.381830\n", " 0.741531 0.279705\n", " 0.949108 0.129485" ] }, "execution_count": 457, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df" ] }, { "cell_type": "code", "execution_count": 458, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Integral exacta 57.0455255009354\n" ] } ], "source": [ "import sympy as sm\n", "import matplotlib.pyplot as plt\n", "from matplotlib.patches import Polygon\n", "from ipywidgets import interactive,interact, HBox, Layout,VBox\n", "import ipywidgets as widgets\n", "from IPython.display import display, clear_output\n", "\n", "\n", "def f(x):\n", " return x*np.sin(x)+x\n", "\n", "a, b = 1,10\n", "\n", "x = sm.symbols('x')\n", "primitiva = sm.integrate(x * sm.sin(x) + x,x)\n", "primitivaN = sm.lambdify(x,primitiva)\n", "integralExacta = primitivaN(b)-primitivaN(a)\n", "print ('Integral exacta ',integralExacta)\n", "\n", "def ErrorRel(aprox):\n", " e=100*abs((aprox-integralExacta)/aprox)\n", " return e" ] }, { "cell_type": "code", "execution_count": 447, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def graficoActualizar(a,b,f,q,k,n):\n", " fig, ax = plt.subplots()\n", " \n", " gap = 0.2 * (b-a)\n", " \n", " fig.set_figheight(10)\n", " fig.set_figwidth(15)\n", " \n", " x = np.linspace(a-gap, b+gap)\n", " y = f(x)\n", " ax.plot(x, y, 'r', linewidth=1) \n", " \n", " fig.text(0.9, 0.1, '$x$')\n", " fig.text(0.12, 0.9, '$y$')\n", "\n", " ax.spines['right'].set_visible(False)\n", " ax.spines['top'].set_visible(False)\n", " ax.xaxis.set_ticks_position('bottom')\n", "\n", " ax.set_xticks((a, b))\n", " ax.set_xticklabels(('$a$', '$b$'))\n", " ax.set_yticks([])\n", "\n", " h = (b-a) / k\n", " \n", " for i in range(0, k):\n", " x0 = a + i*h\n", " x1 = x0 + h\n", " \n", " # Cuadratura\n", " xq = np.linspace(x0, x1)\n", " iq = q(xq,x0,x1,n)\n", " ax.plot(xq, iq, 'k', linewidth=1)\n", "\n", " # Función\n", " iy = f(xq)\n", " ax.fill_between(xq, 0, iq, facecolor='#0099cc', edgecolor ='k')\n", " ax.fill_between(xq, iq, iy, hatch = '//', alpha='0.1') \n", " plt.show()" ] }, { "cell_type": "code", "execution_count": 448, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "ti = {\n", " 0 : np.polynomial.legendre.leggauss(1)[0],\n", " 1 : np.polynomial.legendre.leggauss(2)[0],\n", " 2 : np.polynomial.legendre.leggauss(3)[0],\n", " 3 : np.polynomial.legendre.leggauss(4)[0], \n", " 4 : np.polynomial.legendre.leggauss(5)[0],\n", " 5 : np.polynomial.legendre.leggauss(6)[0],\n", " 6 : np.polynomial.legendre.leggauss(7)[0] \n", "}\n", "\n", "ci = {\n", " 0 : np.polynomial.legendre.leggauss(1)[1],\n", " 1 : np.polynomial.legendre.leggauss(2)[1],\n", " 2 : np.polynomial.legendre.leggauss(3)[1],\n", " 3 : np.polynomial.legendre.leggauss(4)[1], \n", " 4 : np.polynomial.legendre.leggauss(5)[1],\n", " 5 : np.polynomial.legendre.leggauss(6)[1],\n", " 6 : np.polynomial.legendre.leggauss(7)[1] \n", "}\n", "\n", "cNCC = {\n", " 0: [1], #rectángulo apoyado en a\n", " 1: [1/2, 1/2], #trapecio\n", " 2: [1/3, 4/3, 1/3], #simpson\n", " 3: [3/8, 9/8, 9/8, 3/8], #simpson 3/8\n", " 4: [14/45, 64/45, 24/45, 64/45, 14/45], #boole\n", " 5: [5*19/288, 5*75/288, 5*50/288, 5*50/288, 5*75/288, 5*19/288],\n", " 6: [41/140, 216/140, 27/140, 272/140, 27/140, 216/140, 41/140]\n", "}\n", "\n", "cNCA = {\n", " 0: [1], #punto medio\n", " 1: [3/2, 3/2], #trapecio abierto\n", " 2: [8/3, -4/3, 8/3], \n", " 3: [55/24, 5/24, 5/24, 55/24], \n", " 4: [6*11/20, -6*14/20, 6*26/20, -6*14/20, 6*11/20],\n", " 5: [7*611/1440, -7*453/1440, 7*562/1440, 7*562/1440, -7*453/1440, 7*611/1440],\n", " 6: [8*460/945, -8*954/945, 8*2196/945, -8*2459/945, 8*2196/945, -8*954/945, 8*460/945]\n", "}" ] }, { "cell_type": "code", "execution_count": 449, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def NCC(n,a,b,k):\n", " h = (b-a) / k\n", " Q = 0\n", " if (n==0):\n", " for j in range (0, k): \n", " xj = a + j * h\n", " Q += cNCC[n][0]*f(xj)\n", " return Q*h \n", " else: \n", " for j in range (0, k): \n", " q = 0\n", " xj = a + j * h\n", " if (xj==a): \n", " q += cNCC[n][0]*f(xj)\n", " else:\n", " q += 2*cNCC[n][0]*f(xj)\n", " for i in range (1,n):\n", " xi = xj + i * h / n\n", " q += cNCC[n][i]*f(xi)\n", " if (xj + h == b):\n", " q += cNCC[n][n]*f(xj+h)\n", " Q += q\n", " return Q*h/n" ] }, { "cell_type": "code", "execution_count": 450, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def NCA(n,a,b,k):\n", " h = (b-a) / k\n", " Q = 0\n", " if (n==0):\n", " for j in range (0, k): \n", " xj = a + j * h\n", " Q += cNCA[n][0]*f(xj+h/2)\n", " return Q*h \n", " else: \n", " for j in range (0, k): \n", " xj = a + j * h\n", " for i in range (1,n+2):\n", " xi = xj + i * h / (n+2) \n", " Q += cNCA[n][i-1]*f(xi)\n", " return Q*h/(n+2)\n", " \n", "def GL(n,a0,b0,k):\n", " h = (b0-a0)/k\n", " Q = 0\n", " for j in range (0,k):\n", " a = a0 + h*j\n", " b = a + h\n", " q = 0\n", " for i in range (len(ti[n])):\n", " xi = (b+a)/2 + h*ti[n][i]/2\n", " q += ci[n][i] * f(xi)\n", " Q += h*q/2\n", " return (Q)" ] }, { "cell_type": "code", "execution_count": 451, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def lagrange(x,X,Y):\n", " n = len(X)\n", " sum = 0 \n", " for i in range (0,n):\n", " prod = Y[i]\n", " for j in range (0,n):\n", " if (i!=j): prod *= (x-X[j])/(X[i]-X[j])\n", " sum += prod\n", " return sum\n", "\n", "def qgl (x, a, b, n):\n", " if (n==0):\n", " return f((a+b)/2)*(x**0)\n", " else:\n", " X = np.zeros(n+1)\n", " Y = np.zeros(n+1) \n", " for i in range (len(ti[n])):\n", " X[i] = (b+a)/2 + (b-a)*ti[n][i]/2 \n", " Y[i] = f(X[i])\n", " return lagrange(x, X, Y) " ] }, { "cell_type": "code", "execution_count": 452, "metadata": { "scrolled": false, "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "def qncc(x,a,b,n):\n", " if (n==0):\n", " return f(a)*(x**0)\n", " else:\n", " h = (b-a) / n\n", " X = np.zeros(n+1)\n", " Y = np.zeros(n+1)\n", " for i in range (0,n+1):\n", " X[i] = a + i * h\n", " Y[i] = f(X[i]) \n", " return lagrange(x, X, Y)\n", " \n", "def qnca(x,a,b,n):\n", " if (n==0):\n", " return f((a+b)/2)*(x**0)\n", " else:\n", " h = (b-a) / (n+2)\n", " X = np.zeros(n+1)\n", " Y = np.zeros(n+1)\n", " for i in range (1,n+2):\n", " X[i-1] = a + i * h\n", " Y[i-1] = f(X[i-1])\n", " return lagrange(x, X, Y) " ] }, { "cell_type": "code", "execution_count": 459, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "kSlider = widgets.IntSlider(value=1, min=1, max=8, step=1, description=('Intervalos'));\n", "nSlider = widgets.IntSlider(value=1, min=1, max=7, step=1, description=('Puntos (n+1)'));\n", "metodoDropdown = widgets.Dropdown(\n", " options=['Seleccionar Método','Newton-Cotes Cerrado', 'Newton-Cotes Abierto', 'Gauss-Legendre'],\n", " description='Graficar:',\n", " disabled=False\n", ")\n", "\n", "\n", "def GaussInteractivo(n, k, m):\n", " n = n-1\n", " if (m=='Gauss-Legendre'):\n", " graficoActualizar(a,b,f,qgl,k,n)\n", " elif (m=='Newton-Cotes Abierto'):\n", " graficoActualizar(a,b,f,qnca,k,n)\n", " elif (m=='Newton-Cotes Cerrado'):\n", " graficoActualizar(a,b,f,qncc,k,n) \n", " solNCC = NCC(n, a, b, k)\n", " solNCA = NCA(n, a, b, k)\n", " solGL = GL(n, a, b, k)\n", " df = pd.DataFrame({'Integral Exacta': [integralExacta, 0 ],\n", " 'N-C Cerrado': [solNCC, ErrorRel(solNCC)],\n", " 'N-C Abierto': [solNCA, ErrorRel(solNCA)],\n", " 'Gauss-Legendre': [solGL, ErrorRel(solGL)]},\n", " index=['Aproximación', 'Error Porcentual']) \n", " display(df)" ] }, { "cell_type": "code", "execution_count": 460, "metadata": { "scrolled": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "291983aadcf44b45a80a62b852c922b6", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(Dropdown(description='Graficar:', options=('Seleccionar Método', 'Newton-Cotes C…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "widget=interactive(GaussInteractivo, n=nSlider, k=kSlider, m=metodoDropdown)\n", "controls = HBox([metodoDropdown, nSlider, kSlider], layout = Layout(flex_flow='row wrap'))\n", "output = widget.children[-1]\n", "display(VBox([controls, output]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }